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ABSTRACT 



Aims. Despite extensive observations over the last decades, the central questions regarding the power source of the large IR luminosity of Ultra 
Luminous Infra Red Galaxies (ULIRGs), and their evolution, are still not fully answered. In this paper we will focus on massive star formation 
as a central engine and present an evolutionary model for these dust-enshrouded star formation regions. 

■ Methods. An evolutionary model was created using existing star formation and radiative transfer codes (STARBURST99, RADMC and 

RADICAL) as building blocks. The results of the simulations are compared to data from two IRAS catalogs. 

Results. From the simulations it is found that the dust surrounding the starburst region is made up from two components. There is a low optical 
Q ■ depth (r = 0. 1, which corresponds to 0.1 % of the total dust mass), hot (T~400K) non-grey component close to the starburst (scale size lOpc) 
\— i and a large scale, colder grey component (lOOpc, 75K) with a much larger column (r = 10). The simulations also show that starburst galaxies 
C/3 , can be powered by massive star formation. The parameters for this star forming region are difficult to determine, since the IR continuum lu- 

minosity is only sensitive to the total UV input. Therefore, there is a degeneracy between the total starburst mass and the initial mass function 
^> (IMF) slope. A less massive star formation with a shallower IMF will produce the same amount of OB stars and therefore the same amount of 
• »"H . irradiating UV flux. Assuming the stars are formed according to a Salpeter IMF (*P(M) oc M~ 235 ), the star formation region should produce 
' 10 9 M© of stars (either in one instantaneous burst, or in a continuous process) in order to produce enough IR radiation. 
^ , Conclusions. Our models confirm that massive star formation is a valid power source for ULIRGs. In order to remove degeneracies and further 
- - i determine the parameters of the physical environment also IR spectral features and molecular emissions need to be included. 

Key words. Galaxies: starburst - Galaxies: active - Galaxies: nuclei - Infrared: galaxies - Infrared: ISM 



1. Introduction 

In 1983 the Infra-Red Astronomical Satellite (IRAS) surveyed 
96% of the sky in four broad-band filters at 12/im, 25/im, 
60/rni, and lOO/mi. IRAS detected infrared (IR) emission 
from about 25,000 galaxies, primarily from spirals, but also 
from quasars (QSOs), Seyfert galaxies and early type galax- 
ies. Among these galaxies IRAS discovered a new class of 
galaxies that radiate most of their energy in the infrared. The 
most luminous of these infrared galaxies, the (ultra-)luminous 
infrared galaxies [(U)LIRGs], have QSO-like luminosities of 
L > 10 11 Lq (LIRGs) or even L > 10 12 L Q (ULIRGs) 
JGenzel & Cesarsk\l200ol) . 

Despite extensive observations over the last decades, the 
central questions regarding the power source of the large IR 
luminosity of ULIRGs, and their evolution, are still not fully 
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answered. ISanders et alJ (fl988) proposed that most ULIRGs 
are powered by dust-enshrouded QSOs in the late phases of 
a merger. The final state of such a merger would be a large 
elliptical gala xy with a massive q uiescent black hole at its cen- 
ter jKormendv & Sandersl [l992). A significant fraction of the 
ULIRG population seems to confirm this assumption, since 
they exhibit nuclear opti cal emission line s pectra similar to 
those of Seyfert galaxies (Sand ers etalJl988t) . Some also con- 
tain compact central radio sources and highly absorbed, hard 
X-ray sources, all indicative of an active nucleus (AGN). 

On the other hand, the (Far)IR, mm, and radio character- 
istics of ULIRGs are similar to those of starburst galaxies. A 
centrally condensed burst of star formation activity (called a 
starburst, hereafter denoted as SB), for instance fueled by gas 
driven into the center of the potential well of a pair of interact- 
ing galaxies by a b ar instability, prov ides an equally plausible 
power source (e.g. iBlain et alJl2002l) . Observational evidence 
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Fig. 1. Schematic overview of the "physical" environment cre- 
ated in the simulations. In the center star formation is going on, 
which irradiates the dusty surroundings (here represented as a 
toroidal shape). This dust then re-radiates its given energy in 
the IR regime. Also the parameter space is shown. In the star- 
burst, the total stellar mass Msb, the star formation rate (SFR) 
and the IMF-slope a are varied. In the dust the geometry is 
varied by changing the closing angle 77. Also the dust column 
density (expressed as the optical depth r) is varied. A last vari- 
able is the observational inclination S. Also some parameters 
are shown, which initially were not varied: the radius rn and 
the radial scale size dro. More information on the parameter 
space can be found in Sect.|5] 



for the starburst nature of ULIRGs was found with the detec- 
tion of a large number of compac t radio h ypernovae in each of 



the two nuclei of Arp220 by Smit h et alJ dl99{ 

Since both options are observationally supported, it is nat- 
ural to think that there might be an evolutionary relation be- 
t ween the two . Several authors have suggested such schemes. 
In ll988llBaanl explained the evolution of the FIR properties of 
active nuclei with a model which incorporated both a relatively 
rapid decreasing thermal SB component and a slower evolving 
non-thermal component. 

Recently, more papers were published suggesting evo- 
lutio nary schemes wh ich incorporate both SB and AGN. 
iKawakatu et alJ (|2£)06) presented a sample of ULIRGs with 
type I Seyfert nuclei and showed with both observations and 
modeling that this type of galaxies can be powered by both 
a SB and a black hole (BH) with super-Eddington accretion. 
They suggest a scheme in which SB-powered ULIRGs and 
QSOs are two stages in the evolution of the host galaxy. The 
host will start as a SB-driven ULIRG, with only a small BH. 
Over time, the SB fades and the BH will grow into a super 
massive BH (SMBH), which will dominate the energy output 
of the host. By that tim e the host has b ecome a QSO. Simila r 



arguments are mad e by iH opknisetal (20^) and iHol (120 
On the other hand. iKin gl i2005t> . ISilkl J2005I) andlHollEr 



1 

N 
P 
U 
T 

P 
A 
R 
A 
M 
E 
T 
E 
R 
S 



starburst 



SB mass I 
► 

IMF SFR I 



STARBURST99 



r 



stellar radiation field 



geometry 



dust density | 



inclination 



dust 



RADMC 



I 

dust temperature 

* 



RADICAL 



present schemes where the outflow of a super-Eddington ac- 
creting SMBH drives into the surrounding ISM, creating bub- 



Fig. 2. This flowchart shows the computational setup of the 
simulation. First the starburst properties are calculated using 
STARBURST99 (Sect.|2j} and the results of these calculations 
are used as input for the dust calculations, which are performed 
by RADMC and RADICAL (Sects. I2~2landl23l. 



bles in which the gas cools and stars are formed. It is clear that 
detailed modeling is necessary to determine the evolution of 
ULIRGs and their engine(s). 

In this paper, we focus on massive star formation as a cen- 
tral engine. In SB-powered ULIRGs, that are at sufficiently low 
redshift for their internal structure to be resolved, the great 
majority of the IR emission is found to originate from sub- 
kpc dusty regions within merging systems of galaxies (e.g. 
iDownes & Solomonll998l) . These dusty SB galaxies are an im- 
portant class of objects. About 25% of the high-mass star for- 
mation within 10 Mpc distance from us occ urs in just four SB 
galaxies (M82, NGC253, M83, NGC4945; lHeckmarJfl99i) . 
Even though these galaxies create vast amounts of stars, the 
time scale of this formation is short. Near-IR imaging spec- 
troscopy in M82, IC342, and NGC253 indicates that in the 
evolution of these galaxies there are several episodes of star 
formation activity, with timescales of around 10 7 to 10 8 years 
(Genzel & Cesarskv 2000). The relatively low efficiency of the 
energy production of stars (E ~ 10~ 3 M*c 2 ) and the large en- 
ergy output (up to ~ 10 61 ergs for an ultra-luminous starburst 
like Arp220) yield a production 10 8 to 10 10 M© of stars per 
burst. Combined with the short timescales, this leads to star for- 
mation rates (S FRs) ranging from 10 up to 1000 Mq per year 
( Heckman 1998). Numerical simulations confirm this picture. 
Jons son et alJ (2006) show that during a typical merger event, 
there are several short periods of star formation at a very high 
rate. 

The goal of this paper is to present a model for these dust- 
enshrouded star formation regions and to study and explain the 
behavior of the broad band IR continuum properties of star- 
burst ULIRGs during their evolution. The model consists of a 
number of existing codes, which are combined into one. 
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In future work, we will also investigate the spectral fea- 
tures in the IR regime, which will provide more diagnos- 
tics to determine the source of the activity in ULIRGs (e.g. 
Hollenbach & Tielenslll999t iGenzel & CesarskvlkoOOl ISpooiJ 
2003). We also intend to extend the model with a molecular en- 
vironment, in order to further constrai n the physical param eters 
of the cores of active galaxies (Meiierink & Spaans 2005). We 
will also investigate the similarities and differences between 
SB dominated and AGN dominated ULIRGs and will investi- 
gate the possibility of an evolutionary connection between the 
two. 

The structure of this paper is as follows: in Sect. [2] the 
model is discussed. The parameter space and the simulations 
are presented in Sect. [3] and the results of these simulations in 
Sect.|4] In the last section the results are discussed and sugges- 
tions for future work are made. 



2. The model 

The "physical" setup of our model is shown in Fig. \l\ A star 
formation region is placed in the center of our model, sur- 
rounded by a dusty environment, which is responsible for the 
re-emission of the stellar UV photons in the IR. The same 
scheme can be seen in Fig. [21 which shows the computational 
setup of the model. The first step in the simulations is the cal- 
culation of the spectral energy distribution (SED) of the star 
formation region using STARBURST99. The surrounding dust 
is modeled using RADMC and RADICAL, which use the stel- 
lar SED to calculate the dust temperature and the emergent IR 
radiation. The used codes will be briefly discussed and the in- 
terested reader is referred to the listed references for further 
details. 



2.1. STARBURST99 

As mentioned above, the stellar properties of the starburst are 
calculated using STARBURST99 1 (hereafter SB99), which is 
a set of stellar models, optimized to reproduce properties of 
galaxies with active star formation and their evolution through 
time. In this paper only a brief summary of the model is given. 
Mor e information on the data produced by SB99 can be found 
in lLeifherer et al.l d 19991) . Detailed information on t he com pu- 
tational techniques can be fo und in lLeifherer et alJ (l992) and 
iLeitherer & Heckmanl (1 1995b . 

The calculations can roughly be divided into three stages. 
First an ensemble of stars is formed according to the available 
gas mass and the initial mass function (IMF), either in one in- 
stantaneous burst (i.e., the duration of the burst is short com- 
pared to the stellar life time) or continuous with a constant star 
formation rate (SFR). Because these two cases represent the ex- 
tremes of star formation , realistic scenarios will lie somewher e 
in between iLeitherer & Heckmanl 19951 ILeitherer et alJl 992). 

After creating the stellar population, evolutionary models 
are used to calculate the time dependency of the stellar proper- 



ties, such as stellar mass, luminosity, spectra, effective temper- 
ature, and chemical composition. 

In the last stage of the calculation the integrated prop- 
erties of the star formation region are computed by deter- 
mining the stellar number densities, assigning all proper- 
ties to these stars and i ntegrating over the entire population 
(ILeitherer & Heckmanl 19951) . 

2.2. RADMC 

The output SEDs of the SB99 calculations are used to irradiate 
the dusty environment. In order to do this the SB99 spectrum 
must be re-binned to the number of channels the dust simu- 
lation codes can handle. The re-binning of the SB99 output 
is performed in such a way that both the UV/optical part of 
the spectra (which provides the energy for the dust heating) 
and the IR regime (which is the part of the spectrum of in- 
terest) are well covered. Since the SB99 output SED ranges 
from 100 A (extreme UV) up to wavelengths around 10 6 A (100 
//m; FIR) and the number of channels is limited, the resolu- 
tion of the re-binned SED is low. Therefore, the re-binning is 
performed in such a way that the IR part of the SED has a 
somewhat higher resolution, so that the lO/mi, 20/um, 60yum 
and 100/im bands are included, in order to compare the result- 
ing IR data with IRAS observations 2 . Because, in this work, 
we are focusing on defining a physical framework for an evo- 
lutionary model and not on the details, and also taking the 
low spectral resolution and the simple dust model into ac- 
count, we did not compare our results with more recent ISO 
or Spitzer data. Readers who are interested in more detailed 
models which are not taking evolut ion into account are ref - 
ered to, for instan ce, Efstathio u & Rowa n-Robin sonl (^995), 
van Bemmel & Dullem ondf "( l2003l) . iFritz et alJ feOOfj) and ref- 
erences there in. 

The physical dust environment is simulated by creating a 
dust density grid with a central cavity, which is the region sur- 
rounding the central SB where the dust is destroyed and/or 
blown away by the star formation activity. Since we also want 
to explore the influence of the dust geometry on the emergent 
SEDs in our simulations, a 3D grid is needed. This, however, 
would lead to long calculation times. In order to limit this prob- 
lem, the grid is chosen to be axisymmetric. This way the grid is 
reduced to two dimensions again, but still allows for various as- 
trophysically interesting dust distributions, like shells, tori and 
disks. 

RADMC 3 is a 2-D Monte-Carlo radiative transfer code for 
axisymmetric circumstellar disk s and enve lopes. It is based 
on the method of Biorkman & Wood (2001), but with several 
modifications to produce smoother results with fewer photon 
packages. 



1 The code and pre-calculated data sets are publicly available at 
http://www.stsci.edu/science/starburst99/ 



2 In the SB99 output SED, 12 and 25 yum are not available. Instead, 
10 and 20 yum are used, which are good approximations, given the low 
resolution of the data and the simple dust model. 

3 For more information on RADMC, its 
availability and related publications see 
http://www.mpia-hd.mpg.de/homes/dullemon/radtrans/radmc/ 
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Table 1. This table lists the parameters used for the simulations. The top half of the table shows the default simulation and the 
simulations done in the second stage. In this stage the large scale parameters were varied. The bottom half shows the simulation 
from the third stage, where the influence of the second dust component was investigated. All simulations were calculated for in- 
clinations 0, 30, 60 and 90 degrees. The names chosen for the simulations represent the parameter investigated in that simulation. 
Parameters not listed for a simulation were not used in that simulation. 



name 


M SB 


a 


T 


?/ SFR 




dn 77, Ti 


dust 




logM 


M- a 




0-1 M©/yr 


pc pc 


pc 0-1 




DEFAULT 


9 


2.35 


10 


0.3 






grey torus 


LOWMASS 


8 


2.35 


10 


0.3 






grey torus 


HIGHMASS 


10 


2.35 


10 


0.3 






grey torus 


LOWALPHA 


9 


1.3 


10 


0.3 






grey torus 


HIGHALPHA 


9 


3.3 


10 


0.3 






grey torus 


LOWTAU 


9 


2.35 


1 


0.3 






grey torus 


HIGHTAU 


9 


2.35 


100 


0.3 






grey torus 


COVERED 1 


9 


2.35 


10 


0.5 






grey torus 


COVERED2 


9 


2.35 


10 


0.7 






grey torus 


CONTINUOUS 




2.35 


10 


0.3 100 






grey torus 


30pcTORUS 


9 


2.35 


10 


0.3 


30 




grey torus 


50pcTORUS 


9 


2.35 


10 


0.3 


50 




grey torus 


lOpcSHELL 


9 


2.35 


10 


0.3 


10 


10 1 T + T 2 = 


10 " grey torus and grey stellar 


10*10pcSHELL 


9 


2.35 


10 


0.3 


10 


10 1 T+10* 


T2 = 10 * grey torus and grey stellar 


FINAL 


9 


2.35 


10 


0.3 


10 


10 1 0.1 


grey torus and non-grey stellar 


" The total column of both the stellar and the torus dust is scaled to an 


optical depth of 10. See 


Sect.l4.3.2lfor more information. 



b First the density of the stellar dust is multiplied by a factor of 10, then the total column of both the stellar and the torus dust is scaled to an 
optical depth of 10. See Sect. l4.3.2l for more information. 



The method introduced bv lBiorkman & Wood fcOOlh is to 
divide the luminosity of the radiation source (in this case the re- 
binned output spectrum of SB99) into equal-energy, monochro- 
matic "photon packets" that are emitted stochastically by the 
source. The packets are then traced to random interaction lo- 
cations on the dust grid, determined by the optical properties 
and density distribution of the dust. In every iteration step, the 
number of packets absorbed in each of the grid cells is recorded 
and the total amount of injected energy is used to calculate the 
new temperature of that cell. In order to ensure energy con- 
servation and radiative equilibrium, the absorbed energy must 
be re-emitted immediately. The frequencies of the new packets 
are determined by the new grid cell temperature. Local thermal 
equilibrium (LTE) is assumed, so the amount of emitted energy 
equals the absorbed energy. 

Because all the injected packets eventually escape and so 
does all the injected energy, the total energy is automatically 
conserved. Another advantage of this approach is that the simu- 
lation does not have some sort of convergence criterion to meet, 
so that the only source of error in the calculation is the statisti- 
cal error inherent to Mont e Carlo si mulations. 

One problem with the Bior kman & Woodl method is that it 
produces very noisy temperature profiles in regions of low op- 
tical depth, and requires a large number of photons (N ~ 10 7 ) 
for a smooth SED. These disadvantages have been solved in 
RADMC by treating absorption partly as a continuous pro- 
cess (see lLucvlll999l for more information), and using the re- 
sulting smooth temperature profiles with a separate ray-tracing 
code (RADICAL) to produce images and SEDs. These images 



and SEDs have a low noise level even for relat ively few pho- 
ton packages (N ~ 10 5 ). This improved Biorkma n & Wood! 
method works well at all optical depths, but may become slow 
in cases where the optical depth is very large (r v ~ 1000; see 
IPascucci et alJl2004h : this situation is not reached in our simu- 
lations. 

2.3. RADICAL 

RADICAL 4 produces the final spectra, by solving the radiative 
transfer equation on the entire grid, using the dust temperature 
and the SED of the central source to calculate the source func- 
tion. 

First an image is created by integrating the source function 
along rays through the medium. Each ray represents one pixel 
of the image. Spectra can be derived by making images at a 
range of frequencies, and integrating these over an artificial de- 
tector aperture. 

Resolution problems may occur, since the source spans a 
large radial range. The central parts of the image are often much 
brighter than the rest, but cover a much smaller fraction of the 
image. The spectrum may therefore contain significant contri- 
butions of flux from both the central parts and the outer regions 
of the image. Unless the image resolves all spatial scales of 
the object, the spectra produced in such a way are unreliable. 

4 For more information on RADICAL, 

its availability and related publications see 
http://www.mpia-hd.mpg.de/homes/dullemon/radtrans/radical/ 
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Therefore, rather than arranging the pixels over a rectangle, as 
in usual images, they are arranged in concentric rings. The radii 
of these rings are related to the radial grid points of the trans- 
fer calculation. To resolve the central source, some extra rings 
are added. Using these kind of images, all relevant scales are 
resolved, while using only a fairly limited number of pixels 
JPullemond & TuroUj200oT) . 



3. The simulations 

In this section, the simulation setup and the parameter space are 
discussed. As shown in Fig. [2 there are 6 parameters, which 
control the simulation: the combined mass of the stars in the 
starburst Msb, the star formation rate (SFR), the IMF slope a, 
the closing angle of the dust distribution r], the dust column 
expressed in the edge-on V-band optical depth t and the obser- 
vation declination 6. 

This parameter space has been explored, starting with val- 
ues chosen either from literature or at random. In order to verify 
the results of the simulations, a reference data set was made us- 
ing IRAS data of starburst galaxie s taken from the extended 12 
micron galaxy s ample (|R ush et a i1ll993l) and the IRAS Bright 



micron galaxy sample (Kush et al. iy 
Galaxy Sample JCondon etalJl99lh . 



From the simulations a default set is chosen, which fits the 
data best. This set will be referred to as DEFAULT. This set 
was chosen as the starting point for a more systematic sweep 
of the parameter space. Each parameter was both increased and 
decreased with respect to the default value. Exceptions are the 
closing angle rj, the inclination 5 and the SFR. Both permuta- 
tions of rj are higher than the default value. The inclination was 
varied in all simulations: each setup was calculated for 6 = 0, 
30, 60 and 90 degrees. The SFR was only used in few simu- 
lation. In the other simulations, all stars are formed instanta- 
neous. The entire parameter space is listed in Tabled 



3.1. STARBURST parameters 

The SB in the center of the simulation has a large number of 
parameters. Only three are chosen as variables, in order to re- 
duce the amount of data and to avoid unnecessary redundancy. 
These parameters are also expected to have the largest influence 
on the final results. The fixed parameters are listed in Table [2] 

The variable parameters are the total mass of the stars in- 
volved in the starburst (Msb), the star formation rate (SFR) and 
the IMF slope 5 (a). These were chosen because the large lu- 
minosities of ULIRGs can be explained by either a high total 
burst mass, a high SFR or a lower IMF slope exponent. Adding 
more stars or increasing the SFR will scale the entire spectrum, 
whereas decreasing the IMF slope increases the relative num- 
ber of heavy stars, leading to mor e UV radiation a nd therefore 
more heating flux (as suggested in Heckman 1998). 



Table 2. This table lists the parameters which determine the 
simulations, but were not varied. The upper part of the table 
lists the SB99 parameters, the lower part the RADMC parame- 
ters. 



name 


value 




IMF upper mass limit 


100 


M© 


IMF lower mass limit 


1 


M© 


supernova cut-off mass 


8 


M© 


black hole cut-off mass 


120 


M© 


initial time 


0.01 


Myr° 


time step 


0.1 


Myr 


end time 


20 


Myr 


wind model 


theoretical 6 


atmosphere model 


Pauldrach/Hillier* 


metallicity of UV line spectrum 


solar 




Hi 


100 


pc c 


dr 


50 


pc 



" For numerical reasons SB99 cannot start the simulation at t=0. 
h See lLeitherer et all 1 19921) and lLeitherer et all 1 19991) for more in- 
formation. 

c In the third stage, some simulations were made with a different r 



3.2. Dust setup 

The dust distribution is made by creating a polar grid, with a 
density on each grid point. The shape is created by introducing 
a Gaussian distribution in both the r- and ^-direction: 



p(r,ff) - exp 



1 I n/2 -i 
"2\ il 



1 / r - rp 

2 \ dro 



(1) 



The IMF is defined as: Y(M) cc AT". 



where p (r, 0) is the dust density on each grid point, r is the ra- 
dial direction, 9 the angular direction, tq the radius of the dust 
geometry, dro the radial scale size and 77 the closing angle. This 
density distribution is then fixed by scaling the edge-on opti- 
cal depth (obtained by integration along the line of sight) to 
the desired V-band optical depth t. Three examples of result- 
ing density distributions are plotted in Fig. [3] From this two- 
dimensional grid, a three-dimensional grid is made by first mir- 
roring the grid on the r-axis and then rotating this semi-circle 
around the z-axis. 

Two dust components are created: a large scale component 
distributed in a thick disk/torus (here after called torus dust) 
and a small scale shell (called stellar dust). The torus dust dis- 
tribution is created by using an ro of 100 pc and closing angles 
(77) between 0.3 and 0.7. The stellar dust distribution is made 
by taking an ro of 10 pc and an r\ of 1, making it a closed shell. 

Besides the density and geometry, also the optical pa- 
rameters of the dust must be specified. The dust used in 
the s imulations has a grain size distribution n(a) cc cr 3 - 5 
dMathis et al.lll977l) . The torus dust is assumed to be optically 
grey, i.e., radiation of all frequencies is absorbed in the same 
way and re-radiated using a black body profile. The stellar dust 
is much closer to the SB and will be hotter. In order to model 
the expected richer spectral properties of this component, the 
stellar dust is assumed to be composed of graphite, silicates 
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p distribution for simulation: DEFAULT 



p distribution for simulation: COVERED/ 



p distribution for simulotion: 10PCSHELL 






Temperature distribution for simulation: DEFAULT 



Tempe 



distribution for simulation C0VERED2 





Temperature distribution for simulation: 1 0PCSHELL 



Trrin = 35K . 
Tmax = 251K _ 




Fig. 3. These figures show three examples of the dust and temperature distributions resulting from the simulations. The top panels 
show the dust density in grey scale, where dark colors indicate regions with high densities and light colors are low density regions. 
Also some contours are added to enhance the shape of the density profile. The bottom panels show the temperature distribution. 
In these panels, dark colors indicate high temperatures and light colors cooler regions. The minimum and maximum temperature 
is printed in the upper right corner of the panels. Contours are drawn from 50 to 250K in steps of 25K. The data in these plots 
was smoothed over 3 pixels to make the general temperature distribution more clear. The left panels show the default simulation, 
the middle ones a setup with r\ = 0.7 and the right ones a default distribution with an extra 10 pc shell of grey dust. 



and amorphous carbon (e.g. lDraine & LeeTl 984). 

4. Results 

The code produces four broadband IR fluxes as a function of 
time. These fluxes are used to calcu late the IR luminosity (Lir) 
follo wing the definitions g iven in ISanders & Mirabell fl996) 
and lKim & Sanders! Jl998h . 

The simulations were carried out in three stages. First a de- 
fault simulation was made. In the second stage the large scale 
parameters, like the starburst and the torus dust, are explored, 
which influence the long wavelength radiation. The last stage 
focuses on the stellar dust, which affects the short wavelength 
emission. 

4.1. Stage 1 - The default simulation 

The results of the default simulation from the first stage are 
plotted in Fig. [5] The first 4 panels show the results for the four 
IRAS bands. As can be seen the simulation results match the 
reference set quite well. Only the shortest wavelength at 10/mi 
seems to be a poor fit. This is due to the fact that in this stage 
we have not included the stellar dust yet (see Sect. I4.3l for the 



stellar dust). This dust component will radiate mostly at short 
wavelengths. 

The calculated L JR , which is plotted in the fifth panel, also 
covers the range observed in galaxies. At the peak of its activity 
the luminosity of the starburst exceeds 10 12 Lq, making it a 
ULIRG. It remains that luminous for about 3 Myr, after which 
it slowly dims and eventually reaches the LIRG regime. 

In the last three panels of Fig. |5J three IR colors 
(10//m/100jum, 20/im/lOOyi/m and 60/im/lOOyt/m) are plot- 
ted versus the IR luminosity. These color-luminosity dia- 
grams (CLDs) show results similar to the IRAS fluxes. The 
20//m/100jum and 60/ym/100/im results match the observations, 
whereas the lOyum/100/im is a poor fit, due to the absence of the 
stellar dust. 

The default results exhibit a strong dependence on inclina- 
tion. The difference between edge-on and face-on ranges from 
a factor of 3.2 at lOO/zm up to 1000 at 10//m, with an average 
factor of 5.6 in the Lir. This large separation at short wave- 
lengths is caused by the dust distribution. As can be seen in 
Fig. |5J the resulting temperature distribution is not spherical. 
Along the edge-on line of sight the optical depth is 10, whereas 
along the face-on line of sight the star formation region is al- 
most "naked". This causes the areas with lower temperature 
to become more flattened, compared to the inner, high tem- 
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Fig. 4. These panels give an overview of the simulation re- 
sults. The four IRAS bands (in grey: 10/im solid, 20^m dashed, 
60/iin dash-dotted and 100/im dotted) and the IR luminosity 
(black solid) are shown as a function of the various parame- 
ters (top-left: Msb, top-right: a, bottom-left: r, bottom-right: 
77). The points represent the values at 2Myr, which is the time 
where the luminosity of the system is at its maximum. All data 
were scaled to the results of DEFAULT, in order to make an 
easy comparison. Note that the results for the SFR are not in- 
cluded in these plots. Also the observational inclination 6 was 
not included: the data are taken from the face-on results. 



perature areas. When this temperature distribution is viewed 
edge-on, the cooler isothermal surfaces are seen at their max- 
imum surface area and when observed face-on, their surface 
areas decrease rapidly. On the other hand, the inner high tem- 
perature regions are roughly spherical and will therefore look 
the same from all directions. The overall face-on spectrum will 
thus have a larger contribution at short IR wavelengths, as com- 
pared to the edge-on spectrum. This is also illustrated in the 
lower right panel of figure[3] which shows the temperature dis- 
tribution of simulation COVERED2, in which the cooler dust is 
more spherical and the inclination dependence of the lOfim flux 
is an order of magnitude smaller. Another effect of the dust dis- 
tribution is that short wavelength radiation emitted in the center 
will encounter much more dust when it is traveling along the 
edge-on line of sight compared to face-on, making the chance 
of it getting scattered or re-processed much larger. 

The effect of this inclination dependence can also be seen in 
the three CLDs in figure|5] At short wavelengths (10/mi/100//m 
and 20//m/100/mi) the face-on observations trace the data sets, 
whereas the edge-on results cover an entirely different range. 
This effect will be discussed in more detail, when the results of 
the other simulations are presented. 



4.2. Stage 2 - Large scale parameters 

In this section the parameters influencing the star formation and 
the torus dust are explored. Each parameter was varied twice, 
in most cases "in both directions" from the default values, in 
order to uncover systematic effects. The observational inclina- 
tion will not be discussed, since the effects of the inclination 
have already been discussed above. The figures which show 
the most important effects have been included in this article. An 
overview of the results is shown in Fig. [4] The rest of the fig- 
ures displaying all the results, including the four IRAS fluxes, 
of the individual simulations and comparisons between them 
are collected in Appendices A and B respectively. 

4.2.1 . Starburst mass 

The effect of the starburst mass on the SB99 spectra is quite 
simple. Since this code works with a distribution of stars, inte- 
grated to a certain mass, the spectra will scale proportional to 
the total stellar mass. The influence on the dust is less straight- 
forward, however. When the energy input provided by the star- 
burst increases, the overall dust temperature will be higher and 
the flux at short wavelengths will increase faster relative to the 
long wavelengths. Also the IR luminosity is expected to in- 
crease with the stellar mass, since the energy input increases. 

The simulations confirm these expectations. Figure |4] 
which gives an overview of the results of the simulations, 
clearly shows that the short wavelengths (grey solid and 
dashed) increase much faster than the long wavelengths (grey 
dash-dotted and dotted) and that Lir (black) increases as the 
mass increases. 

The color-luminosity plots in Fig.[6]seem to reveal less triv- 
ial information. Assuming all the other parameters are correct, 
a minimal SB mass of 10 9 Mq is needed in order to create an 
IR luminosity higher than 10 12 Lq. Also all three simulations 
seem to lie on the same "track" in the CLDs, and starburst mass 
and age have the same effect on the results. In other words, a 
young, but relatively light starburst will, in IR colors, have the 
same appearance as a massive, older starburst. 

4.2.2. IMF slope 

When the IMF slope decreases, the number of massive stars in 
the starburst will increase. These high mass stars are the main 
source of UV radiation, and so the dust temperature will in- 
crease with a shallower IMF slope. This means that varying the 
IMF slope would yield results similar to increasing the starburst 
mass. 

Figure |4] shows this is indeed the case. Note that at first 
sight the effect seems to be reversed, but the IMF is defined as 
^(M) oc M~ a , and the axis is reversed. Interestingly, the effect 
is somewhat more complicated than for Msb- All fluxes scale 
proportionally to Msb, but the effect is stronger when going 
from high a to the default value, then when lowering it even 
further. This is caused by the fact that the number of OB stars 
does not increase linearly with a. 

Figure B.2 shows that although LOWALPHA starts with 
higher values for the IR luminosity and the colors, after about 
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10 Myr the low IMF case cannot be distinguished from the de- 
fault simulation any more. This is due to the short life time 
of the high mass stars. Therefore, the same conclusion can be 
drawn for the IMF slope as for the starburst mass: the data ei- 
ther represent young starbursts with a normal IMF, or an older 
population of stars with a shallow IMF slope. 

4.2.3. Dust column density 

Several effects have been identified when changing the total 
amount of dust for a given energy input. 

First, when the dust column is lower, the dust will be 
heated more evenly, since the radiation can travel through the 
dust easier. This will result in a smaller temperature gradient, 
with a higher minimum temperature. The spectrum will be- 
come flatter, and have less contribution from 10;um and 20fim. 
Conversely, large dust columns lead to larger dust temperature 
gradients and steeper spectra. This effect is shown in Figure 0] 
where the variation of t has effects similar to Msb and a. 

Second, due to the larger temperature gradient, the inclina- 
tion dependence increases for larger t. In Figure this effect 
can be seen in the first panel, where the edge-on and face-on 
results are almost the same for LOWTAU, whereas the results 
for HIGHTAU differ by up to an order of magnitude. 

A third effect is the increase in the total IR luminosity Lir, 
simply because there is more radiating dust and therefore more 
UV photons can be converted into IR. 

All these effects can be found in the color-luminosity dia- 
grams. When going from LOWTAU to DEFAULT especially 
Lir increases, resulting in a shift to the right. But when t in- 
creases even more, the color changes rapidly, and the track 
moves to the upper right corner of the plot. Again, due to the 
inclination effects, the edge-on results do not match the refer- 
ence sets. The edge-on results of LOWTAU are better, because 
they suffer less from inclination effects, but they give worse fits 
than the face-on DEFAULT results. 

Like in the previous simulations, there is a degeneracy be- 
tween t and time. Young starbursts surrounded by moderate 
column density dust have the same IR properties as older, more 
enshrouded starbursts. 

4.2.4. Closing angle 

The effect of the closing angle is not as profound as that of the 
other parameters. The only major effect is that the inclination 
dependence effects are reduced, which is caused by a smaller 
density gradient in the 9 direction (also see Fig. [3}- The extreme 
case would be a simulation with T] = 1, which is a shell of dust 
with no inclination dependence at all. 

Although the inclination effects are less in simulations 
with a higher closing angle, the edge-on tracks in the color- 
luminosity tracks of Fig.[8]still do not fit the reference data. 

The shape of the spectrum also changes when the closing 
angle changes, but this is due to the fact that the total dust col- 
umn changes, since the edge-on optical depth is kept constant. 



4.2.5. Star formation rate 

So far, all simulations were made under the assumption that 
the stars were formed instantaneously, or, in other words, on 
a timescale much shorter than the evolutionary timescale of 
the stars (< 1 Myr). This is based on the assumption that there 
is a rapid build-up of gas in the starburst environment, which 
at some point becomes critically dense and "ignites". It also 
assumes that the stellar winds and supernovae do not provide 
strong feedback, so that all gas can be converted into stars on 
a short time scale. However, the assumption that the starburst 
environment is less violent and more like normal star formation 
regions, with a steady inflow of gas, is also valid and must be 
explored as well. 

In order to address this issue, simulations were done with 
continuous star formation scenarios with different SFRs and 
durations as presented in Figs. A. 10, B.5 and B. 10. It was found 
that a SFR of several tens of M0 yr~' is required to reproduce 
the fluxes found in observational data. Simulations with a du- 
ration of 100 Myr show that after about 10-20 Myr there is 
no further evolution in the IR fluxes and colors. However, our 
data cannot exclude longer burst durations, but there are con- 
siderations that would make such scenarios unlikely. First of 
all, long periods of star formation at high rates would require 
vast amounts of gas. Secondly, these high star formation rates 
also produce a lot of stellar feedback (stellar winds, SNe) that 
makes it increasingly difficult to supply gas to form new stars. 

Because of these reasons, we consider the case of a high 
SFR and a short duration (100 M0 yr _I and 10 Myr, leading 
to the formation of 10 9 M0 of stars, and making it compara- 
ble to the other simulations). The results of this simulation are 
compared to DEFAULT in Fig. B.5. This figure shows some re- 
markable differences between the two situations. Whereas the 
instantaneous burst starts at a maximum and decays rapidly af- 
ter about 3 Myr, the continuous case starts at much lower val- 
ues and takes quite some time to rise. After about 5 Myr the 
two scenarios cross and the continuous case starts to dominate. 
After about 10 Myr it reaches its maximum, but it does not 
decay, since stars keep forming at a constant rate. 

The significant difference between the continuous and the 
instantaneous star formation is that objects with instantaneous 
star formation are only luminous enough in the first 5 to 10 
Myr, which decreases their detection rate. The continuous star- 
bursts remain luminous as long as there is inflow of gas to feed 
the star formation. 

4.3. Stage 3 - Stellar dust 

The torus dust alone is not able to explain all the IR contin- 
uum properties and therefore in the second set of simulations 
addresses the properties of the stellar dust component. The lack 
of emission at short wavelengths can have two reasons. A very 
simple explanation is that the geometry we are currently using 
does not provide enough hot dust and more dust must be added 
in the center of the dust distribution. Another possibility is that 
the optical properties of the large-scale dust are inappropriate 
and that the stellar dust modelling requires more realistic opti- 
cal properties. 
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4.3.1 . Smaller radius 

One of the ways to create more hot dust is to decrease the size 
of the total dust distribution. This way all the dust will be closer 
to the central source and therefore the overall temperature will 
be higher. In order to investigate this possibility, two simula- 
tions were carried out, with ro = 50 pc and rn = 30 pc. The 
rest of the parameters were the same as for DEFAULT. As can 
be seen in Fig. B.6, decreasing ro does improve the short wave- 
length results. However, it also affects the longer wavelengths 
and although the fit in the lO/mi/lOO/vm CLD is better, the slope 
is not correct. A positive effect is that is introduces scatter in 
the vertical direction of the CLD, which could not be explained 
well by the previous simulations. Overall one can conclude that 
although changing the radius of the dust geometry does im- 
prove the results, it does not provide the hot dust contribution 
that we are looking for. 

4.3.2. Grey stellar dust 

Another possible solution is to add additional stellar dust close 
to the starburst. Maybe our assumption that the dust close to 
the central region is completely destroyed or blown away is not 
correct. Therefore, we introduce an extra (Gaussian) dust shell 
with a scale size of 10 pc (see also Fig. 01. Two simulations are 
performed, one where the total column of both the stellar and 
the torus dust is scaled to an optical depth of 10 and one where 
the dust density of the stellar component is first multiplied by 
a factor of 10 and then the total column is scaled to an opti- 
cal depth of 10. Again, the rest of the parameters are the same 
as for DEFAULT. Fig. B.7 shows that introducing this second 
component does increase the IQ/im flux, without affecting the 
60yt/m and 100//m fluxes. However, the 20jum flux is still influ- 
enced, and the 10/im flux cannot be increased further, without 
compromising the 20 fim results. These simulations also intro- 
duce scatter perpendicular to the evolutionary tracks. 

Both attempts to solve the lO/vm problem by introducing 
extra hot dust, do not provide the results that are desired. This 
suggests that the optical properties of the dust must be modi- 
fied. 

4.3.3. Non-grey stellar dust 

Since the stellar dust component must have non-grey opti- 
cal properties, it it will be composed fr om graphite, silicate s 
and amorphous carbon (as described in Drain e & Lea ll984). 
In these computations, the SB99 output spectrum is first pro- 
cessed by the non-grey stellar dust and the emergent radiation 
is used to irradiate the (grey) torus dust distribution. 

As can be seen in Fig.[9](see also Fig. A. 1 1 and Fig. B.8 for 
the full results and Fig.[l0|and Fig. A.12 for SEDs), the addi- 
tion of non-grey stellar dust finally produces the right amount 
of 10/mi radiation, without influencing the flux at longer wave- 
lengths. In addition, although this second dust component has 
a large influence on the lO/zm radiation, only a very low optical 
depth is needed, which corresponds to a small amount of dust 
(see Fig. B.9). The final results are created using an V-band op- 
tical depth for the second component (T2) of only 0. 1 . In terms 
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Fig. 10. These panels show the evolution of the emergent spec- 
trum of our final simulation. The top panel shows four spec- 
tra as a function of time ranging from t=0 (black) to t=15Myr 
(light grey) in steps of 5Myr. The bottom panel shows the same 
spectra, but this time divided by the spectrum of t=0. All spec- 
tra were taken from the face-on results. 



of mass this means that the amount of dust in the second com- 
ponent is about 0.1% of the total dus t mass. Observati onally, 
similar values are found. For instance, Klaa s et alJ(ll997l) found 
that in order to properly fit the SEDs of Arp 244, NGC 6240, 
and Arp 220, they need to add about 0.1% of hot dust to the 
bulk of the (colder) dust. 

Physically, the stellar dust component most likely repre- 
sents the left-overs of the dust, which was shocked by stellar 
ejecta and blown away by the stellar radiation. This dust is 
much closer to the stars than the torus dust, is hotter and experi- 
ences disruptive events like sputtering and shattering. A simple 
grey dust model does not give the right results, since for these 
smaller grain sizes the wavelength depend ence of the absorp- 
tion c ross section is important (e.g. S poonlEo03T : IPeeters et alJ 
2004J. The torus dust is the ring or shell of dust surrounding the 
starburst nucleus. It has higher column densities, larger grain 
sizes and is further away, giving rise to lower dust tempera- 
tures. Therefore, the grey approximation is valid for this dust 
component. 

5. Conclusions 

Our simulation model is able to reproduce the IR continuum 
properties of starburst galaxies, by using two dust components: 
a large scale component containing the bulk of the mass and a 
less massive component close to the starburst. The long wave- 
length radiation is influenced by macro-physics like the star 
formation activity and the large scale dust distribution. The 
short wavelength emission at 10 /im comes from hot dust and is 
influenced by the stellar dust and micro-physics like the optical 
properties and optical depth of the dust. 
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Not all parameters have a profound effect on the results. 
The stellar dust optical depth fo) and the parameters control- 
ling the star formation (the total stellar mass Msb and the IMF 
slope a) have a large effect on the final results, whereas the 
influence of the large scale dust geometry (r and rf) on the final 
IR properties is smaller. 

Star formation region 

Because of their large influence on the results, the star for- 
mation parameters could be determined well. Increasing the 
stellar mass (Msb) has two effects. First of all, the shape of 
the spectrum changes, the short wavelength fluxes increase 
as compared to the long wavelengths. A second effect is an 
increase in Lir, which changes a factor of 10 when changing 
the mass by a factor of 10. The effects of varying the IMF 
slope (a) are similar to those of the stellar mass, but are 
less pronounced when decreasing a than when increasing it. 
Making the IMF slope more shallow increases the luminosity 
by a factor of 4.0, whereas steepening the slope decreases the 
luminosity by a factor of 7.9. There is a degeneracy between 
Msb and a. A less massive starburst with a shallower IMF 
will produce roughly the same amount of OB stars as a more 
massive starburst with a steeper IMF and therefore the same 
amount of massive stars and therefore the same irradiating UV 
flux. Assuming the stars are formed according to a Salpeter 
IMF (*P(M) oc M~ 2 35 ), the star formation region should 
produce 10 9 Mq of stars (either in one instantaneous burst 
or in a continuous process) in order to produce enough IR 
radiation. 

Stellar dust 

For the stellar dust component, the grey approximation of 
the optical dust properties is not valid. A more realistic dust 
model, including graphite, silicates and amorphous carbon, is 
necessary to produce the right results. In addition, the V-band 
optical depth (tz) is found to be an important factor. Using a 
T2 of 10 (or 1) overproduces the amount of 10/im radiation by 
almost a factor of 4 (or 2.5). We require an optical depth of 0. 1 
to get values that agree with observations. 

Torus dust 

The influence of the large scale dust geometry on the final 
IR properties was far less and therefore these parameters 
could not be determined to great precision. Increasing the dust 
column density (t) has two effects. Like for Msb, the short 
wavelengths are enhanced compared to the long wavelengths 
and the IR luminosity increases. These effects are, however, 
much smaller: the luminosity changes by a factor of 2.0 when 
going from t=1 to 10 and by a factor of 3.2 when increasing 
it further to 100. A second effect is that the inclination 
dependence increases with increasing optical depth. The 
difference in the edge-on and face-on values of Lir changes 
by a factor of 1.6 for t = 1 to a factor of 16 for r = 100. The 
best fit to the data was obtained with r = 10, but the other 
simulations were also reasonable. The closing angle did not 
seem to have an optimum value at all. Varying it only affects 
the inclination dependence of the results. The variation of Lir 
for a flat disk-like structure (77 = 0.3) is about a factor of 6.3, 



compared to a factor of 2.0 for a more shell-like dust geometry 
with a closing angle of 0.7. Even though the inclination effects 
are reduced for higher values, the edge-on results still do not 
match the data. 

Considering the values determined for the parameters 
investigated, it seems that, although observationally starburst 
galaxies appear to have a very violent nature, the star formation 
environment does not need to be as "exotic" as one might 
expect. The only exceptional parameter needed to explain the 
high IR output of starburst galaxies is the large amount of 
massive stars (high Msb or SFR), but there is no need for an 
adjusted IMF to increase the number of heavy stars. The torus 
dust surrounding the stars is no exception to this trend. Both 
the size (100 pc) and the V-band optical depth of the torus 
dust are moderate (10, comparable to values found for photon 
dominated regions). Furthermore, most of the dust has a low 
temperature (~ 75K at lOOpc from the center), while only a 
very small amount of hot dust (~ 400K at lOpc) is needed. 

In all simulations, including the final ones, the data were 
only fit well by the face-on evolutionary tracks. All the com- 
pletely edge-on results were a poor fit. This effect is stronger in 
the short wavelength results than in the long wavelength result. 
This indicates that inclination dependence is mostly caused by 
the obscuration of the hot dust in the center by the outer dust 
distribution. This can have two implications. On the one hand, 
it could be that the inclination dependence resulting from our 
model is too large, and that extra parameters are needed to ad- 
dress this problem. A likely parameter is the dumpiness of the 
dust. A clumpy medium with the same average density (i.e. the 
same mass) would hav e a lower apparent optical d epth then 
a smo oth medium (e.g. iNatta & Panagial[l984l IConwav et al.l 
2005), which makes it easier for the short wavelength radiation 
to travel in the edge-on direction. On the other hand, if the pre- 
dictions of our model are correct, the implication is that there is 
a significant number of starburst ULIRGs, which are currently 
not classified as such based on their IRAS colors. 

A specific shortcoming of our models is that, although the 
IR properties are well explained, almost all parameters move 
the evolutionary tracks more or less along the same line and in 
the direction of the evolution. The result is that not all parame- 
ters of the physical environment in a given starburst galaxy can 
be uniquely inferred from observations, using this model. The 
UV input can be inferred from the total IR output, but this does 
not constrain the IMF or the SFR. Similarly, the optical proper- 
ties of the dust can be determined, but the geometry is hard to 
infer. To address these problems, more information than just the 
IR continuum is needed and therefore we intend to extend the 
current model. First of all, the IR part of the code will be mod- 
ified to include specific spectral characteristics (e.g. PAHs and 
high ionization lines). Also the molecular environment that sur- 
rounds the current dust region will be added to the model. The 
chemistry of such a region will give a better handle on the radi- 
ation fiddj_as_w£ll_asJh£^nsiJdes and temperatures of the gas 
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Also other wavelengths will be studied, since optical and UV 
data will give more information about the star forming re- 
gion, whereas (sub)millimeter and radio observations will re- 
veal more about the outer dust regions and the molecular envi- 
ronment. 
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Fig. 5. These panels show the infrared results of the simulation DEFAULT. The top 4 panels show the 10, 20, 60 and lOOfim 
flux as a function of time. The fifth panel shows the calculated infrared luminosity (Lir) through time. The last three panels give 
a more observational interpretation of the data. In these panels three IR colors (10pm/100/im, 20/mi/100jum and 60yum/100/im) 
are plotted versus the IR luminosity. To indicate the evolution direction, the start of the track is marked with a star. All axes are 
logarithmic except the time axes. The observational inclinations are indicated by the line style: from face-on in solid, through 
dashed and dash-dotted to edge-on indicated with dotted lines. The two reference data sets are added with markers. In the first 5 
panels, the erro r markers represe nt the mean of these data sets and a 1 cr deviation. The gre y marker is data fro m the 12 micron 
galaxy sample (Rush et al. 1993) and the black one from the IRAS Bright Galaxy Sample (iCondon et a l. 1991). In the last three 
panels, the data from the 12 micron galaxy sample is indicated by circles and the data from the IRAS Bright Galaxy Sample by 
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Fig. 6. These panels show the influence of the starburst mass (Msb) in the infrared results. The two deviations are compared to 
the default value. The simulation results for the lower deviation are indicated by dotted lines, the higher deviation by dashed 
lines and the default set is drawn in dash-dotted lines. The face-on observations are marked with a star, the edge-on results with 
a circle. These marks also indicate the start of the evolutionary track, to indicate the evolution direction. The panels are the same 
as the bottom four panels of Fig. [5] the first panel shows the IR luminosity as function of time and the other three panels show 
three IR colors (lO/im/lOOpm, 20/im/100^/m and 60pm/100/im) versus the IR luminosity. As in Fig.|5]the two reference data 
sets are added with markers. In the first panel, the error mark ers represent the m ean of these data sets and a 1 cr deviation. The 
grey ma rker indicates data fr om the 12 micron galaxy sample jRush et alJl993l) and the black one from the IRAS Bright Galaxy 
Sample lCondonetalll99ll> . In the last three panels, the data from the 12 micron galaxy sample is indicated by circles and the 
data from the IRAS Bright Galaxy Sample by squares. The black dots connected by a solid line show the reference data binned 
per 10 points. 
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Fig. 7. These panels show the influence of the dust optical depth (r) on the infrared results. The color coding is the same as in 
Fig.0 
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Fig. 8. These panels show the influence of the closing angle (77) on the infrared results. The color coding is the same as in Fig. [6] 
Note: for this simulation both deviations are higher than the default value. 
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Fig. 9. These panels show the influence of adding the stellar dust component on the infrared results. The color coding is the same 
as in Fig. [6] Note: to show that the whole reference set is traced in the CLDs, 20 Myr is shown instead of 10. 



